Two-Species Reaction-Diffusion System with Equal Diffusion Constants: 
Anomalous Density Decay at Large Times 
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We study a two-species reaction-difFusion model where A + A — > 0, A + -B — > and B + B — > 0, 
with annihilation rates Xo,So > Ao and Aq, respectively. The initial particle configuration is taken to 
be randomly mixed with mean densities nA(0) > ns(0), and with the two species A and B diffusing 
with the same diffusion constant. A field-theoretic renormalization group analysis suggests that, 
contrary to expectation, the large-time density of the minority species decays at the same rate as 
the majority when d < 2. Monte Carlo data supports the field theory prediction in d = 1, while 
in d = 2 the logarithmically slow convergence to the large-time asymptotics makes a numerical test 
difficult. 
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Fluctuation effects in low-dimensional reaction- 
diffusion systems has received considerable attention 
lately [Q . When the dimension d of the system is suffi- 
ciently small, below a critical dimension dc, spatial fluc- 
tuations in the particle concentration dramatically influ- 
ence the large-time behavior of the system. Most strik- 
ingly, the large-time scaling exponents and, sometimes, 
amplitudes of various observables become universal and 
independent of the microscopic description of the system. 
Indeed, the quest for possible new universality classes 
has intensified in recent years, with studies of a vari- 
ety of models exhibiting fluctuation-dominated kinetics 
at different levels of complexity. Theoretical advances, 
together with a growing number of experimentally ac- 
cessible realizations of low-dimensional reaction-diffusion 
models - like reaction kinetics of excitons in polymer 
chains or surface deposition of certain proteins ||l| - have 
further accelerated these studies. 

We here consider a generalization of the well-known 
single-species A + A ^^ (I) and two-species A + B ^ (it 
annihilation models f^ by coupling them together and 
allowing all particles to react: 



A + A 



A + B 



So 



B + B 



\'n 



(1) 



The A and B particles are assumed to perform a random 
walk on a d-dimensional lattice with diffusion constants 
Da and Db respectively. When two particles meet on 
the same site they annihilate with probabilities given by 
the reaction rates Aq, Sq, and Ag in (|l|). Among a num- 
ber of possible applications, the model may be used to 
describe steric reaction-diffusion processes ||3|, with the 
two species representing the projections of the chemically 
"active spots" of the reactants onto some fixed reference 
direction Q. 

The question we wish to address is how the average 



particle densities nA(t) and nsit) decay in time when 
there is an initial imbalance nA{0) > ns(0) in the pop- 
ulation of the two species. (In an application to a steric 
reaction-diffusion process, this would correspond to a 
bias in the initial orientations of the reactants, triggered 
by some external field Q.) When 5q > Ao,Ao we ex- 
pect the minority species of _B-particles to die out be- 
fore the majority A. This has been confirmed by a 
field-theoretic renormalization group analysis for the case 
when Da 7^ Db f^]: Independent of dimensionality, the 
minority species was shown to vanish before the extinc- 
tion of the majority. Here we shall focus on the case 
when Ag = Aq and with the diffusion constants tuned to 
the same value, Da = Db = D, as would be the proper 
choice when applying the model to a steric reaction- 
diffusion process [Bj. Since the effective time-dependent 
reaction rates are controlled by the diffusion constants 
B, it is interesting to explore whether the case of equal 
diffusion constants could lead to new effects. 

To set the stage, let us briefly study the classical rate 
equations implied by (1): 



da 2 

-di = -(^°" 



Soab), ^ = -{Xob^ + doba), (2) 
dt 



with a, b the mean-field densities of species A, B, and 
with initial conditions a(0) = 71^(0) and b{0) — ns(0). 

The qualitative behavior of (^) can be studied by a 
mapping to a Poincare sphere |q] with z — 1/a and 
u = b/a: 



dz ^ r ,,-, du 
-^Xo + 6ou{t), - 



(Ao - SoHl - u)/z , (3) 



where u is assumed to be in the interval [0,1]. With 
•^0 < So, (S) implies that for u(0) < 1 and for large times 
l/u{t) diverges as 



u{t) b{t) 



oo, 



(4) 



when a(0) > 6(0). Thus, with Ao and (5o fixed, the minor- 
ity species is killed off early, as expected. Note, however, 
that (0) suggests that u{t) could saturate to a constant 
if the effective diffusion- controlled reaction rates replac- 
ing Aq and (5o in the presence of fluctuations approach 
each other asymptotically. To explore this possibility it 
is convenient to pass to a field-theoretic formulation of 
the problem. 

For this purpose, let us consider the master equation 
of the system: 



dt 



P{c, t) = J2 Rc'^cP{c', t)^Yl Rc^c'P{c, t). (5) 



Here P{c,t) denotes the probability of a configuration c 
at time t, and Rc'^c is the transition rate from state c' 
into c, determined by D^Sq and Aq. We take P(c, 0) as 
a Poisson distribution, with averages denoted by n^fO) 
and 715(0) for A and B particles respectively. Eq. (||) 
can be translated into a Schrodinger-type equation with 
a second-quantized Hamiltonian H , which in the contin- 
uum limit turns into a field theory with action S: 

S= fd-'x fdtiY, Mdt-'^^)<f>o.+ J2 ^"-s^ 

•' •' a=A,B a.,P=A,B 

x^,,{\ + 4>fi/2)(j)^(j)p]~ I d'^x Y, "a(O)0a(a^,O), (6) 

•^ a=A,B 

with rescaled time variable, Dt -^ t. The scalar fields 
(j)A and 4>B describe the A and B particles with initial 
densities ^^(O), a — A,B, and with the coupling con- 
stants Xaa — Xbb = Ao and Xab = Aba — 5q encoding 
the reaction rates. The action in (||) faithfully mimics 
the master equation in (0), up to terms irrelevant under 
renormalization. Note that the rate constants become 
dimensionless when d = 2 which suggests that dc — 2, a 
result confirmed by standard power counting S. Since 
we are interested in effects from fluctuations, we focus on 
d < 2 in what follows. 

When d < 2 the particle densities can be calculated 
perturbatively from (g) as na{t) =< (pa >s using an 
e— expansion, with the rate constants Aq and Sq replaced 



by renormalized rates gain ) 



^ao/(l 



'ao/5* 



with a = X and (5, respectively, and with k the renormal- 
ization scale. Divergent logarithmic terms in the expan- 
sion are conveniently grouped together using "running 
coupling constants" 



5A(i) = 



9 



1 + (t/n) 



e/2-' 



gsit) = 



1 + (t/ti)-^/^ 



(7) 



where g* = r(e/2) ^(Stt)'^'^ is the renormalization group 
(RG) fixed point to which g\(t) and gsit) flow under a 



change of time scale Q . The running coupling constants 
in (|7|) are obtained from the Callan-Symanzik equation 
which expresses the independence of the densities Ua (t) 
on the choice of k: 

na{t;g\iK~^),gs{K^^);nA,o,nB,o;K) = {K^ty^/"^ x 
xnc(K''^;gx{t),gs{t);nAfi{t),nBfl{t)]n) , a ^ A, B (8) 

with fiafiit) = {n^tY/'^nafi. 

Employing the mean-fleld densities from ( B) in (@) 
(fluctuation-improved mean-field treatment [p[|10[), we 
obtain 



L{t) -fu* = (r2-i)/(fi + i) , 



(9) 



where Vi^ = 1 + 4^*uo/(l - mf and T = exp(-87r(3i^ - 
J-)), 1 for d = 2, 1. Thus, at this level of analysis the 
minority species at large times decays at the same rate 
as the majority [u* ^ 0) . Note also that (0) implies that 
u{d = l,t) > u{d = 2,t) 'AS t -^ CX3, suggesting that 
the survival rate of the minority particles increases with 
reduced dimensionality. 

The fluctuation-improved mean-field analysis (where 
all tree-level diagrams are summed) takes into account 
effects from fluctuations via the renormalized coupling 
constants. However, in order to fully account for fluctu- 
ations one has to explicitly consider the loop diagrams. 
For this purpose it is convenient to expand in the asymp- 
totically small parameter r](t) = g\{t) — gs{t). To one- 
loop level and to first order in 77 we thus obtain 



na{t)^^laN{t){l + Aa{t)) , 



A,B. 



(10) 



u*/(l + u*) and 
00, and we hence 



Here ^ia = 1/(1 + "*), A*s = 
A A (i) , Ab (i) ex Tj{t) In i'*/^ ^ as i 
recover the result in (|9|) at very large times. The common 
factor N{t) in (nO) coincides with the density of the single 
species A + A — > model, with N{t) = Int/Snt md=2 
and N{t) = l/VS^t'^/'^ in d = 1 0. The approach 
of the total density nA{t) -\- nsit) to that of the single- 
species model reflects the fact that the equal values of 
the diffusion-controlled rates in (|^) at large times makes 
the two species indistinguishable, effectively leaving us 
with a single-species system. We should here stress that 
the values of the amplitudes in (ll^) have been obtained 
by ignoring all irrelevant terms (in RG sense) that could 
be added to the action in (|6|). Since we do not expect 
the amplitudes to be universal, such terms - if included 
- could shift the value of u* in (^). The scaling form 
of nA{t) and nsit) on the other hand is expected to be 
universal (considering the generic behavior of this class 
of problems 0), and should be insensitive to any left-out 
irrelevant terms. 

To summarize: Our RG analysis to one-loop level sug- 
gests that the minority species at large times decays at 



the same rate as the majority when d < 2. This result 
- if robust against contributions from higher-order terms 
in the loop expansion - is quite remarkable and very dif- 
ferent from the conventional picture implied by classical 
rate theory. Unfortunately, due to the difficulty of keep- 
ing track of possible "correction-to-scaling" terms g , the 
present model (in contrast to the simpler A + A ^ (d 
model) is intractable to standard procedures for assessing 
higher-loop effects O . Therefore, to check the robust- 
ness of the result in (1^ we have resorted to a Monte 
Carlo simulation of the master equation (pi) (which in 
addition provides information on the transient kinetics). 

To perform the simulation we have used an improved 
minimal process algorithm, [l2[ : At each time step a lat- 
tice site, say i, is picked at random (out of L possible 

sites), and a table of statistical weights Wa for the pos- 
sible processes tTq, a = 1,...,M, at this site is con- 
structed, with 7rj\/ a "null process" . The weights are de- 
fined as cumulative probabilities: Wa = X^fci ^0 /Qi 

with Q — X]a=i Ra y where Ra denotes the rate for ttq. 
Q is the total rate which is kept the same for all sites by 
properly adjusting the rates KfJ j — I, . . . ,L for the null 
processes. After selecting a random number < r < 1 
that process tTq is carried out for which w^_i < r < Wa 
{w_[ = 0), and time is incremented hy t -^ t + 1/LQ. 

Fig. 1 shows results for the density ratios nA{t)/nB{t) 
from simulations in 1, 2 and 3 dimensions. Initially the 
system is homogeneous with randomly mixed A- and B- 
particles, with n^(0) > ^^(O). At early times a typical 
minority particle is then surrounded by several A's, and 
is therefore likely to be killed off. As seen in Fig. Il|, this is 
indeed the case: in all dimensions the density ratio grows 
initially. However, for sufficiently large rates there is a 
recovery effect: After the initial growth the density ra- 
tio drops — the minority species finds a way to recover. 
(Here and in the following "recovery" is understood in 
relative terms. The number of minority particles still 
decays in time but the ratio of majority to minority den- 
sities decreases in favor of the minority). As seen in Fig. 
1 the large-time behavior of the curves depend on the 
dimensionality as well as on the values of the reaction 
rates. For 6 = 3Ao = 50 the d = 3 curve takes off after 
the recovery phase, in agreement with the classical mean 
field result (Q). This is in striking contrast to the d = 1 
curve which saturates to a plateau, as predicted by the 
one-loop RG result in ([10|). The eventual fate of the d = 2 
curve is more difficult to foretell since the convergence to 
the asymptotic scaling regime in d = 2 is logarithmically 
slow nw. Unfortunately, with available techniques it is 
hard to push the d > 2 simulations reported in Fig. 1 
much beyond logio(i) ~ 0.5 where the number of parti- 
cles get too small to obtain a reliable statistics. Turning 
to the 6 = 3Ao = 10 runs, the d = I curve again saturates 
to a plateau, while the d — 2 and d = 3 curves take off 
even before a peak has developed. For this case it may 



be harder to argue that the d — 2 curve may eventually 
also saturate, although it cannot be excluded considering 
the slow convergence. Note that the RG prediction of an 
enhanced survival rate in low dimensions is confirmed for 
this choice of reaction rates as well: the d—1 curve falls 
below that for d — 2. The runs with very large reaction 
rates, 5 = 3Ao = lO'^, produce curves which collapse onto 
a single curve, spuriously suggesting a plateau in all di- 
mensions: In d = 3, however, large reaction rates impede 
the crossover to the mean field asymptotics, and we ex- 
pect the d = 3 curve to take off at a time not accessible 
with our simulation. On the other hand, we do expect the 
d = 1 curve to stay flat (considering the generic behavior 
for smaller rates), while the asymptotics of the d — 2 
curve is again hard to predict due to the logarithmically 
slow convergence to the scaling region. 

The build-up of a density-ratio peak with a subsequent 
"recovery phase" in the early stage of the kinetics (for 
sufficiently large rates) is made transparent by solving 
the master equation (0) with the diffusion constant D 
set to zero (i.e. with an ensemble of decoupled single site 
problems p^). The hatched curves in Fig. 1^ represent 
the corresponding density ratios for very large and small 
reaction rates. Initially the dominant A + B — > reac- 
tions on the different sites cause a fast rise of nA{t)/nB{t) 
(peak formation). After some time this process gets ex- 
hausted with single B-particles left at a finite fraction of 
the sites. The A + A — > reactions on the remaining 
sites now become dominant, causing a depletion of A- 
particles (recovery phase for B particles). As this process 
in turn gets exhausted one is left with a constant ratio 
nA{t)/nB{t) (plateau). As seen in Fig. ^, for small rates 
the early-time decoupled-sites kinetics gets perturbed by 
diffusion of particles from neighboring sites, which leads 
to an enhanced mixing of A and B particles and hence a 
steeper rise of the density-ratio curve. Still, in d = 1, the 
curve eventually saturates to a plateau. This plateau, 
however, is of a very different origin from that of the 
decoupled-sites problem with no diffusion. This is strik- 
ingly shown in Fig. where we have plotted the individ- 
ual densities in d = 1 for large and small rates, with the 
corresponding density ratio curves inserted at the top. 
In both cases the individual densities are seen to decay 
after the plateau has formed. Note in particular that the 
decoupled-sites plateau for large rates (with no decay of 
the individual densities) persists only in the time inter- 
val —2.5 <logio(t) < —1. At later times effects from 
diffusion sets in, without, however, visibly distorting the 
plateau. Note also that although the asymptotic curves 
come very close to each other, they do not collapse onto 
a single curve, in agreement with the one-loop RG re- 
sult: The total density nA{t) + nB{t) in (pii| ) is universal, 
whereas the individual components are not. 

It is instructive to look at a Monte Carlo "cartoon" of 
the kinetics in d = 1, as plotted in Fig. pi. The "recovery 
phase" (for B particles) is visible as dense red ropes of 
A particles terminating in black (empty sites), with the 



remaining particles diffusing as single threads colored red 
(A particles) or blue (B particles). Remarkably, at large 
times these threads self-organize spatially in such a way 
as to keep the density ratio at a constant value. By a 
comparison to the A + B ^ model EW one observes 
that the tendency for domain formation in the present 
model is much weaker: A and B particles appear to be 
rather well mixed asymptotically. This makes the anoma- 
lous density decay even more intriguing, suggesting that 
the average A-A and A-B distances are correlated asymp- 
totically so as to precisely compensate for the smaller 
number of A-B pairs (with endangered B particles) as 
compared to A-A pairs. 

To conclude, we have studied a two-species reaction- 
diffusion model where the particle densities exhibit an 
anomalous behavior at large times: The spatial self- 
ordering of the two species allows the minority population 
to persevere, with the asymptotic density ratio (of ma- 
jority to minority species) kept at a constant value. This 
apparent "non-Darwinian" kinetics is entirely driven by 
fluctuations in the diffusion-controlled regime, typical for 
a system of reduced dimensionality. A perturbative RG 
analysis suggests that the same remarkable behavior is 
present also in two dimensions, although this is still to 
be verified by independent methods. Details and exten- 
sions will be published elsewhere [|ll| . 
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FIG. 1. Density ratio nA{t)/nB{t) obtained from a numer- 
ical simulation of the master equation (fel) in d = 1, 2 and 3 on 
a 10^, 10^ X 10^ and 10^ x 10^ x 10^ lattice, and numerical solu- 
tion of (H) in the limit of D = ( "decoupled-sites problem"). 
The fixed parameter values are: D = ls~^ , nA,s(0) = 2, 1 
particles/site. The reaction rates Ao and So are given in units 
of the diffusion constant D. 
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FIG. 2. Plot of individual densities (bottom) together with 
density ratios nA{t)/nB{t) (top) in d = 1 obtained from sim- 
ulations. 




FIG. 3. "Cartoon" of the d = 1 kinetics. Time runs hori- 
zontally, with 300 lattice sites represented along the vertical 
axis. At the top the corresponding density ratio is plotted 
for comparison. The reaction rates are Xo,So = 10,30. Red 
(blue) color denotes that a majority of particles belong to 
the A (B) species. White indicates the same number of A's 
and B's on a lattice site, with intermediate cases covered by 
intermediate colors. Black color denotes an empty site. 



